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ABSTRACT 

We model the constraints set on the evolution of dusty starburst galaxies by 
the current deep extragalactic surveys performed in the far-infrared with Spitzer, 
and at radio wavelengths with the VLA. Our models fit the number counts in 
all the available spectral bands well, and also provide a reasonably close match 
to the redshift distribution of the Spitzer detections. We find: 1.) dusty star- 
burst galaxies with infrared burst phases triggered by galactic interactions at 
redshift z ~ 1 — 2 are good candidates to fit the Spitzer results at 24/xm, 70/im 
and 160/zm, assuming plausible strengths for the PAH features for the infrared 
luminous sources. An Arp220-like spectral energy distribution (SED) for Ultra- 
luminous Infrared Galaxies (ULIGs) of L ir > 10 12 L G and one like that of M82 for 
Luminous Infrared Galaxies (LIGs) of L iT ~ 10 11-12 L & give a successful fit to the 
Spitzer 24/im and ISOCAM 15/zm number counts at flux levels of S v < 1 mJy; 
2.) the strong evolution of the number density of the ULIGs from redshift z ~ 
to ~ 1 predicted by our models is consistent with the current deep 1.4 GHz radio 
surveys and accounts for the upturn in the 1.4 GHz differential counts at the sub- 
mJy flux level; and 3.) comparing with number counts at near infrared bands, 
as well as the background measurements using DIRBE and 2MASS, shows that 
only a fraction of the stellar mass in the Universe is included in our models of 
dusty starburst mergers at z ~ 1 — 2. 



Subject headings: evolution-galaxies-interaction-galaxies-starburst-galaxies-Seyfert 
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1. Introduction 

During the last few years, remarkable progress has been made in probing the distant 
Universe. This progress is driven in part by the availability of powerful new instruments at 
various wavelengths spanning the X-ray, near-, mid- and far-infrared, and submillimeter as 
well as out to the radio, which have complemented the traditional studies of galaxy evolution 
at z > 1 based on UV/optical observations. 

A crucial advance is the discovery by far-infrared and submillimetre observations of a 
significant population of dusty starburst galaxies at high redshift with infrared luminosities 
L ir > 10 12 L©, which are ULIGs with central activities (AGN or starburst). These objects are 
heavily hidden by dust extinction and are triggered either by violent mergers between gas rich 
objects or by encounters of galaxies (Smail et al. 1997; Hughes et al. 1998; Blain et al. 1999; 
Eales et al. 2000; Holland et al. 1999; Puget et al. 1999; Dole et al. 2000; Hughes 2000; Sanders 2000; 
Sajina et al. 2003; Xu et al. 2004). The formation rates for massive stars (M > 5 M Q ) are 
approximately 200 — 900 M Q ?/r _1 , assuming the bolometric luminosities of ULIGs are dom- 
inated by star formation activities (Condon 1992; Scoville et al. 1997; Ivison et al. 1998; 
Ivison et al. 2000). 

Another important clue to galaxy evolution is the population of optically faint, near- 
infrared bright Extremely Red Objects (EROs) often found as counterparts of sources iden- 
tified in many other wavebands (Richards 1999; Smail et al. 1999; Alexander et al. 2001; 
Chary k Elbaz 2001; Cowie et al. 2001; Lutz et al. 2001; Pierre et al. 2001; Smith et al. 2001; 
Alexander et al. 2002; Franceschini et al. 2002; Ivison et al. 2002; Mainieri et al. 2002; Page et al. 2002; 
Stevens et al. 2003). The classifications of this population are divided between dusty star- 
bursts and more evolved objects. Morphological studies by Hubble Space Telescope (HST) 
and ground-based near-infrared (NIR) imaging show that active dusty EROs are mostly dis- 
turbed systems, indicating mergers or galactic interactions could be the driving mechanisms 
for their behavior (Smail et al. 2002b; Gilbank et al. 2003; Moustakas et al. 2004). The de- 
reddened star formation rates in these objects are in a range of 20 ~ 200M Q yr~ 1 , with 
far- infrared luminosities generally below 10 12 L & . In these cases, near-infrared deep surveys, 
with detection sensitivity at least three magnitudes deeper than that of far-infrared and 
submillimeter surveys, can provide a possibility to unveil the high-z dusty starburst galaxies 
in a way complementary to the surveys at FIR/submillimeter, which are only sensitive to 
the very bright infrared sources (Kennicutt 1998; Cimatti et al. 1999; Cimatti et al. 2002; 
Smail et al. 2002; Bouwens et al. 2005). 

NASA's new infrared facility Spitzer, which brings significant improvements in IR sen- 
sitivity and survey efficiency over the previous deep infrared cosmological surveys, has pro- 
duced deep new multi-wavelength number counts from 24/im to 160/im that provides a 
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unique opportunity for galaxy evolution studies. Because of their sensitivity to the promi- 
nent PAH spectral features (Polycyclic Aromatic Hydrocarbons), Spitzer as well as the 
Infrared Space Observatory (ISO) 15/im deep cosmological surveys can provide an unbi- 
ased view on the evolution of both ULIGs and LIGs at the redshift range < z < 2.5 
(Egami et al. 2004; Le Floc'h e tal. 2004; Marleau et al. 2004). Although the initial phe- 
nomenological models fitted the unified ISO /Spitzer number counts well, they have been 
shown not to fit the source redshift distribution (Perez-Gonzalez et al. 2005). This paper 
uses models based on a binary aggregation description of galaxy evolution to fit the new 
data and to provide a more complete test of our understanding of infrared galaxy evolution. 

In Section 2 of the paper, we will review our previous work. Section 3 discusses the 
main features of the binary aggregation dynamics and the numerical simulation approach 
we use, which is centered on an infrared burst phase triggered by gas-rich mergers. A brief 
discussion of the construction of template spectra of dusty starburst galaxies is also given in 
this section. The model fitting to the recent Spitzer results, as well as other available multi- 
band deep surveys spanning from near-infrared to radio wavelengths, and the estimation 
of the integrated background intensities are all presented in Section 4. In section 5, we 
summarize our conclusions. The cosmological parameters H = 70km/s/Mpc, Q = 0.3 and 
A = 0.7 are adopted throughout. 



2. Summary of Previous Work 

Theoretical modelling of the evolution of dusty starburst galaxies constrained by both 
the observed Cosmic Infrared Background and the number counts in far-infrared and submil- 
limeter deep surveys by IRAS, ISO, and SCUBA indicates a dramatic increase in the number 
density of LIGs and ULIGs towards redshift z ~ 1 (Dole et al. 2000; Wang & Biermann 2000; 
Chary & Elbaz 2001; Franceschini et al. 2001; Chapman et al. 2002; Totani & Takeuchi 2002; 
Wang 2002; Lagache et al. 2003). Wang & Biermann (2000) and Wang (2002) discussed the 
effects of galaxy mergers on the strong evolution seen in these deep surveys in the context of 
a binary aggregation evolutionary scheme. The model requires that the infrared behavior be 
luminosity- and redshift-dependent with a significant population of dusty starburst galaxies 
or AGNs from gas rich mergers at redshift z ~ 1 — 2, and less luminous gas-poor objects 
at lower z. The infrared emission of more massive merging systems is taken to be enhanced 
to a higher level and to fade away faster than less massive systems within the merger time 
scale of a given epoch. This hypothesis is based on the observation that ULIGs are usually 
more than a factor of 20 brighter than normal starburst galaxies and it is consistent with 
the recent proposal for a "downsizing formation scenario" . 
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The mechanism driving the strong evolution of infrared luminous sources is still unclear. 
It could be, for example, an increment of major/ minor merger rates, an upturn in the 
probability of minor tidal interactions, or an enrichment of cold gas toward high redshift. 
Although visual morphological classifications of merger pairs is still very challenging, work by 
Bell et al. (2005) and Shi et al. (2005) shows that a decrease of the major-merger rate since 
redshift z ~ 1 is probably not the dominant cause for the rapid decline in the cosmic star 
formation rate. This finding supports the prediction by Wang & Biermann (2000), that sub- 
mJy sharp upturn in the far-infrared number counts diagram could not be interpreted solely 
by a merger rate decrease with cosmic time, but required a luminosity-dependent infrared 
burst phase from gas rich mergers at high redshift and gas poor mergers at low redshift. We 
hope that future sophisticated morphological studies from deep imaging as well as CO mass 
determination from ALMA will allow such issues to be addressed more completely. 

Although the details of the infrared burst phase are still unclear, it is believed to be 
related to a stage of the merger process when the dust mass and temperature are both 
dramatically increased (Kleinmann & Keel 1987; Taniguchi & Ohyama 1998). Numerical 
simulations on the evolution of dusty starburst galaxies by Bekki & Shioya (2001) show that 
there is very strong photometric evolution during the merger process of two gas-rich disks, 
and a dramatic change of the SED around a cosmic time scale T ~ 1.3 Gyr, when the 
two merging disks become very close and suffer from violent relaxation. The star formation 
activities at this moment may reach a maximal rate of ~ 4OOM yr _1 . The infrared flux in 
this case increases by an order of magnitude, especially in the far-infrared (60 \xm ~ 170 \xm) 
in the emitting frame. 

However, such dramatic changes only occur for major mergers. Episodic starburst 
phases are indicated by recent observations (Shi et al. 2005), and they may tend to smooth 
out the extreme maxima predicted in the models. Our model is not restricted to major 
mergers. However, for simplicity, we assume in the current work that galactic interactions 
would trigger only a single infrared burst within a merger time scale and simulate a range of 
such interactions with a luminosity-dependent infrared burst phase in addition to the mass- 
light scaling relation of normal starburst galaxies. We believe that episodic bursts would not 
change the main results of this work, but will relax the requirements for a high fraction of 
mergers at high redshift. Thus, the current work should be indicative of the overall pattern 
of infrared galaxy evolution regardless of the detailed time dependence of the star formation. 

The redshift distribution of ULIGs in the model of Wang & Biermann (2000) and Wang 
(2002) shows a remarkable increase of the number density from the local Universe up to 
redshift z ~ 1 and a mild decrease towards still higher redshifts. We have found a similar 
result based on the study of the cosmic star formation history using X-ray deep surveys 
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(Wang et al. 2003). This paper examines whether the evolution of the ULIGs observed 
in the IRAS and ISO deep far-infrared surveys is still compatible with recent observations 
such as: 1.) in the 1.4 GHz radio band, which is itself a star formation indicator, and is 
an important way to study the evolution of starburst galaxies; and 2.) 24/im, 70/im and 
160/im observations by Spitzer, which are particularly useful for studying the evolution of the 
starburst population in the redshift range 1 < z < 2. Another important rationale for this 
work is to try to learn more about the relevant physical processes. Although near-infrared 
number counts, unlike mid- or far-infrared counts, would not provide strong constraints on 
the star formation activities in our modelling, we include here a comparison of our modelling 
with the results of near-infrared deep surveys, to get an estimate of the fraction of stellar 
mass contained in these starburst merging galaxies to the total stellar mass budget in the 
Universe. 

With model constraints ranging from near-infrared out to deep radio surveys, we have 
constructed an overall model for both ULIGs and LIGs that evolves in a continuous way, 
based on a merger triggered starburst scenario. This approach differs from other models 
that assume ULIGs evolve separately from LIGs and with an extremely strong evolutionary 
rate. 



3. Model Description and numerical simulation 

3.1. Key Numerical Relations 

We describe the modelling algorithm briefly in this section. Readers are referred to Wang 
& Biermann (2000) and references therein for the details of the binary aggregation dynamics 
and the numerical techniques. An aggregation phenomenon based on the Smoluchowski 
(1916) equation is adopted in the model to trace the evolution of the mass distribution 
function with cosmic time, which is in the form of: 



<9N(m,t) 
dt 



-y rm poo 

- / dm K(m , m — m ,t)N(m ,t)N{m — m ,t) — N(m,t) / dm K{m,m ,t)N(m ,t) (1) 
2 Jo Jo 



N(m, t) is the mass distribution function in the "comoving" frame. The kernel K(m, m\t) = 
n g (t) < TV [t] > reflects the interaction rate for each pair of masses (m,m), and depends 
on the mechanism and environment of such aggregations. n g is the density of galaxies; V(t) 
is the relative velocity of the interacting pairs, and £ is the cross section. In aggregation 
dynamics, the interaction kernel K = n g < SU > is the key point for driving the whole 
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evolutionary process, which depends strongly on the environmental structures. However, the 
determination of the interaction kernel suffers from various uncertainties and complexities. 
We shall adopt in this model a simplified formula for the aggregation kernel K(m, m ,t) with 
separated time evolution and cross section terms as below: 



K(m,m ,t) oct k (m 2/3 + m' 2/3 ) 1 + a(m 2/3 + m' Z/a ) . (2) 

where t k reflects the evolution of the aggregation rate with cosmic time. The free parameter 
k depends on the specific structures and interacting environments. A range of values is 
discussed by Cavaliere & Menci (1997) for groups, clusters and large scale structures. The 
second free parameter is a, which describes the relative importance of two kinds of encounters 
(geometric collisions and focused resonant interactions) in the cross section. The aggregation 
time scale r is inversely proportional to the interaction kernel K(m, m \t) i.e., K(m, m ,t) oc 

A Monte Carlo inverse-cascading process is used to simulate the binary aggregation 
of galaxies described by Eq. (1). This equation has been solved in a number of cases 
(Safronov 1963; Trubnikov 1971; Silk & White 1978), which show that the dynamical pro- 
cess does not strongly depend on the initial conditions, i.e. the mass spectrum is indepen- 
dent of initial details with self-similar evolution. For this study, we adopted a 5-function at 
M = 2.5 • 10 9 M Q as an initial mass spectrum (so dwarfs are included in the new calculation) 
and start the simulation from redshift z = 15. In the previous studies, we had set the initial 
mass to M — 2.5 • 10 10 M Q . Other model parameters of the dynamical processes are discussed 
in detail by Wang & Biermann (2000) and are still valid for the current modelling. 

Up to now, the mass distribution function N(m, t) could be derived from the Smolu- 
chowski equation using the Monte Carlo simulation. However, to compare the modelling 
with various observational constraints, such as luminosity functions of galaxies in general 
or with the morphologies, source counts, redshift distributions and background intensities 
from multi-band surveys, we need to understand the conversion of the mass distribution 
function N(m, t) to the observable luminosity function N(L, t), as well as the modelling of a 
luminosity-dependent infrared burst phase from gas-rich mergers, all of which are described 
in detail by Wang & Biermann (1998, 2000) and Wang (1999, 2002). We summarize several 
main points as follows. 

A simple prescription of the mass-light ratio for the faint blue galaxies is given by Cava- 
liere & Menci (1997) as L/L* = (M/M*) v (L* is the local standard characteristic luminosity 
with the corresponding mass M#). Here, rj = 4/3 if the cross section is purely geometrical; 
this value is consistent with the observational results by Kormendy (1990) . oc f(z, A , O ) 
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can be used to describe redshift dimming or luminosity evolution. Simplifying the color and 
K- corrections, we get roughly that L B oc ^% M v . We assume in our model an infrared to op- 
tical color ratio L ^ m oc (i) ~ 1.5 is adopted in the calculation). This value is consistent 
with both the numerical calculations of the photometric and spectroscopic evolution of dusty 
starburst galaxies, and with the current understanding of the ULIGs, which are believed to 
be the extremely luminous infrared burst phase due to starburst merger events where the 
far-infrared luminosity, L ir , is enhanced both by accumulation of dust and by an increase in 
the dust temperature. This burst phase can enhance the infrared luminosity by a factor of 
about 20 over that of normal starburst galaxies (Kleinmann & Keel 1987; Silva et al. 1998; 
Taniguchi & Ohyama 1998; Bekki & Shioya 2001). We adopt L m ^ m oc f( z ) AT**, 
f(z) oc (1 + z) 13 in the modelling. A value of (3 ~ 5 after a transition redshift z ~ 1 
gives the best fit, which successfully interprets the sharp upturn at the faint end of far- 
infrared number counts. In addition, choosing z ~ 1 as a transition redshift, with gas rich 
mergers at z > 1 and gas poor mergers at z < 1, is also consistent with the cosmic time 
scale at z ~ 1, which is about 3 x 10 9 years, approximately the time scale for disk evolution. 
Galaxies may become gas poor during this stage. The scaling factor of the mass-light ratio 
is normalized by the local luminosity function from the IRAS survey. 

Except for understanding the physical basis for the two major free parameters in the 
interaction kernel (k ~ —1.35 and a ~ 1), which we adjust in the simulation to give a good 
fit to the current observations, the uncertainties in the relative emission over the radiation 
spectrum across cosmic time as galaxies evolve, including progressive evolution and episodic 
starburst activities triggered by galactic interactions etc., are now the major challenges to 
galaxy evolution studies. With only number counts, redshift distributions, and integrated 
background levels, which are rather loose constraints, we simply assume a luminosity depen- 
dent infrared burst phase after a transition redshift, to mimic the change of the SED shapes 
(especially in the far-infrared) for galaxies with different luminosities. We simulate such 
an effect by adding a scaling factor (^ L )~ < ' to the luminosity evolution term defined in the 
last paragraph (f(z) oc (1 + z) 13 ), where L c is the minimal luminosity bin of the calculated 
luminosity function. The best fit result gives ( ~ 0.9 at far infrared wavelengths. 

This luminosity dependent infrared burst reflects a physical reality, that the infrared 
luminous galaxies at the bright tail of the luminosity function become gas poor faster than 
the less luminous ones and fade away quickly within the merger time scale of that epoch. 
This behavior is consistent with the current concept of a "downsizing formation scenario" . In 
addition, it indicates that less luminous infrared sources, such as LIGs, could have infrared 
burst phases persisting longer than those of ULIGs. In fact, we can describe such an idea 
with an infrared burst "duty circle," f on . This term is based on expressing the duration 
of the infrared burst phase in terms of the merger time scale, t s b/t mer ger, where t s b is the 
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starburst time scale and t mer ger the merging time scale. Learning from observations, we 
have t merger oc (1 + z)~ m and t s b oc Therefore, we can derive an empirical relation 

ir 

fon oc (1 + z) m . From the discussion above, we now understand that this luminosity 

ir 

dependent infrared burst phase is actually implicitly related to the starburst time scale. 

3.2. Model Characteristics 

While the basic concepts of our model are clearly plausible, the details of the necessary 
assumptions are somewhat arbitrary. However, it seems that they are very important for 
fitting the sharp upturn at faint flux levels in the current deep far-infrared surveys. For a 
detailed comparision with the color behavior vs. flux densitiy, we shall in the future include 
a real physical model including both proper chemical evolution and dust emission properties. 

In the current calculations, starburst galaxies, dust shrouded AGNs and "non-evolving" 
galaxies are treated separately as three major classes of infrared emitting sources. The "non- 
evolving" galaxies represent the big spiral and elliptical galaxies, which were already in place 
about 8 Gyr ago and are considered to have a fixed luminosity function with cosmic time 
(Lilly et al. 1998; Schade et al. 1999). We know from recent work by Hammer et al. (2005) 
that some intermediate-mass spiral galaxies may have experienced starburst phases due to 
galactic interactions, which are counted within the starburst population triggered by galactic 
mergers/interactions in our calculation. Although the treatment of the three populations 
separately in the modelling is a bit arbitrary, this simplification should not affect the current 
results significantly. 

Unlike other models that assume ULIGs evolve separately from LIGs and with an ex- 
tremely strong evolutionary rate, we prefer to think that merger-triggered infrared luminous 
sources including both ULIGs and LIGs exhibit a continuous evolution, where the infrared 
emission of ULIGs may change more rapidly than that of less luminous sources (i.e. LIGs 
may have a longer infrared burst duration than that of ULIGs). We simulate such a con- 
tinuous range of evolution by assuming a luminosity dependent infrared burst phase based 
on the consideration of different time scales for the starburst galaxies of different infrared 
luminosities to consume the cool gas. The AGN component is included as an additional 
model constraint, since the dust enshrouded geometry could result in significant radiation 
from them at infrared wavelengths. We assume AGN activities are triggered by galactic 
mergers or interactions based on the aggregation dynamics, similar to the starburst pop- 
ulation. The statistics of the starburst and AGN populations are normalized by the ob- 
served local luminosity function of starburst galaxies at 60 \im from Saunders (1990) and 
that of Seyferts from Rush et al. (1993). The abundance of dust-shrouded AGNs is very 
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uncertain, especially for those dust enshrouded QSOs at high redshift. For the models, 
we set the abundances of the dust-shrouded AGNs relative to the total AGN population 
to be 50% and 80% at local and high redshift, respectively, based on the statistics from 
Hubble Space Telescope imaging survey of nearby AGNs and the cosmic X-ray background 
(Malkan et al. 1998; Gilli et al. 1999). 

3.3. Template SEDs 

To construct a template spectrum for the high redshift dusty starburst sources around 
redshift z ~ 1, we utilize the publicly available program GRASIL introduced by Granato 
and Silva to produce a library of SEDs. We replace the PAH features with a template for 
PAH emission based on rest-frame ISOCAM/CVF 5 — 16.3/im spectra of Arp220 and M82 
(Weedman et al. 2004). We selected two template SEDs respectively for ULIGs and LIGs, 
which are adequate to fit the present deep surveys, especially the 24/im number counts by 
Spitzer and the ISOCAM lbfxm results. A comparison of the template SED for ULIGs and 
LIGs with the measurements of Arp220 and M82 is given in Fig. 1, which shows that the 
dust re-radiation peaks around 50 — 60 fim, consistent with a dust temperature of T ~ 55 K 
and a dust absorption coefficient K\ oc A _/3 with the coefficient index j3 ~ 1.2. A third 
template, the SED of M51, is assumed for the spiral population. 

An important change has been made in our current modelling, compared with our earlier 
work (Wang 2002). Previously, only a simple template SED (without specific consideration 
of PAH features) was assigned for the whole infrared luminous starburst population, rather 
than the separate SEDs for LIGs, and ULIGs. To give a good fit to the new Spitzer results, 
it requires adjustment of the merger rate especially at higher redshifts and of the free param- 
eters in the prescription of the infrared burst phase. This probably contributes significantly 
to the better fits we achieve for the redshift distributions of the infrared sources compared 
with previous efforts that have assumed a single SED (e.g., Lagache et al. 2004). 

We choose the template spectrum for the obscured AGNs at low redshift to be that of 
NGC 1068, and the unobscured ones that of the mean SED of a Seyfert I. The early phases 
of these AGNs are assumed to show typical spectra such as the dust shrouded F10214+4724, 
and a phase poor in cold gas like the Cloverleaf quasar. The templates of all these spectra 
were modelled by Rowan- Robinson (1992), Rowan- Robinson et al. (1993) and Granato et 
al. (1994, 1996, 1997). 
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4. Results and discussion 

4.1. Mid- and Far-Infrared Number Count Fitting 

We first show in Fig. 2 and Fig. 3 the comparison of our model results with the number 
counts from recent Spitzer MIPS deep surveys at 24/im, 70/im and 160/xm. The green dashed 
line is used to denote the population of starburst galaxies with an infrared burst phase due 
to galactic interactions. The black dot-dotted line corresponds to the non-evolving galaxies, 
and the blue dash-dotted line is for Seyferts/AGNs. The sum of the three populations is 
represented by the red line. A similar set of lines will be used in all the following plots. 

The left panel of Fig. 2 is the model fitting of the differential source counts in Spitzer 
24/im deep surveys, which include ~ 5 x 10 4 sources to an 80% completeness of ~ 80fiJy. 
Spitzer's increased sensitivity and efficiency in large area coverage over previous infrared facil- 
ities, coupled with the encompassment of PAH emission in Spitzer's 24/xm band, dramatically 
improve the quality and statistics of number counts in the mid-infrared (Papovich et al. 2004). 
The right panel of Fig. 2 shows the calculated contribution to the differential number counts 
at 24/im from ULIGs and LIGs separately, as well as the sum of the two populations. As 
shown in Fig. 3, we find that the faint end upturn of the number counts at 70/im and 160/xm 
can be interpreted successfully by the population of starburst mergers from our modelling 
with infrared luminosities L 1T > 10 12 L Q and Arp220-like SEDs. 

Fig. 4 and Fig. 5 present the model fitting of the number counts from ISOCAM 
15/xm deep survey, IRAS 60/xm, ELAIS 90/im final analysis and FIRBACK 1 70/im surveys. 
Fig. 4 (Right) to Fig. 6 (Left) show that the number counts for a reliable subset of the 
detected sources at FIR/Submm wavelengths could be sufficiently accounted for by the 
infrared burst phase when a population of ULIGs with Li r > 10 12 L Q might be produced by 
merger-triggered starburst/AGN activities at z ~ 1 (Kawara et al. 1998; Elbaz et al. 1999; 
Efstathiou et al. 2000; Dole et al. 2001). However, we also learned from the simulation that 
the FIR slope of the adopted SEDs would significantly affect the FIR/Submm number count 
fitting at faint flux levels, while the free parameter ( influences mostly the bright part of 
the source count diagram. Similar changes as for the Spitzer data were made in the other 
far infrared number count fittings, basically for the same reasons. For example, the peaks 
of the calculated differential number counts at IRAS 60/im and ISO 90/im are shifted to a 
fainter flux (from ~ 50 mJy to ~ 10 mJy) (the exact observational values are still open for 
debate due to the detection limits of current far- infrared surveys). 

The simulation shows that the extremely strong evolution at sub-mJy of the Ihjim and 
24/im number counts is mainly dominated by LIGs with enhanced PAH emission (see Fig. 
2 (Right)), while ULIGs contribute mostly to the bright part of the peak at 0.3 mJy < S v < 
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1 mJy. The bright end of the number counts diagram is mostly contributed from 

"non-evolving" spiral population and dusty AGNs. We also see that the PAH spectrum 
and strong silicate absorption of the ULIGs affect significantly the contribution of the dusty 
starburst galaxies to the Spitzer 24/im and ISO 15/im number count diagrams in the flux 
range 1 mJy < S v < 10 mJy. This is one of the reasons why the starburst population is 
almost negligible above a few mJy in the newly calculated 15/im number counts, while it 
was dominant up to ~ 10 mJy in the previous work (Wang 2002). Of course, there are a 
few other modifications to better fit the recent Spitzer 24/xm number counts: 1.) we include 
the dwarf population with initial mass spectrum M = 2.5 x 1O 9 M ; 2.) we increase the 
scaling factor of the interaction kernel K(m,m ,t) in equation (2) to shorten the merger 
time scale of dwarfs or minor mergers at high redshift to a reasonable range and enlarge the 
absolute value of the free parameter k to make the exponential damping faster; and 3.) the 
free parameters in the prescription of infrared burst duration ((3, () are adjusted to change 
the relative infrared burst duration of LIGs and ULIGs. 



4.2. Number Counts at Other Wavelengths 

The model fitting of the differential radio counts at 1.4 GHz is shown in Fig. 6 (Right). 
We found that the well known 1.4 GHz source excess observed at sub-mJy levels can be 
interpreted satisfactorily by the population of starburst galaxies with an infrared burst phase 
at z ~ 1. These objects are LIGs of L ir ~ 10 n ~ 12 L Q with enhanced PAH emission. We 
believe that the 1.4 GHz source excess and the strong evolution detected by ISOCAM 15/im 
and Spitzer 24/im deep surveys at flux S u < 1 mJy are due to the appearance of the same 
population. 

Fig. 7 (Left) shows the model fitting of the cumulative number counts of galaxies at 
1.6/j.m brighter than a given flux threshold F M . The data are NICMOS H-band observations 
covering 1/8 of the area of the Hubble Deep Field North that reach a flux level down to 
~ nJy (about three magnitudes deeper than other near- or mid-infrared deep surveys). The 
number count fittings of the K-band and ISOCAM 6.7/zm deep surveys are presented in Fig. 
7 (Right) and Fig. 8 respectively. We know from these results that the population of dusty 
starburst merging galaxies in the model only contributes to the bright end of the number 
counts at K-band and 6.7/zm. There is a significant deficiency of sources from our modelling 
at the faint flux levels (< 0.01 mJy), which is clearly shown in Fig. 7 (Left). This indicates 
that deep near-infrared observations reveal a population of very faint and high redshift 
galaxies that are not detected in current mid-infrared surveys. We estimate roughly that 
30% to 50% of the stellar mass in our Universe could be contained in the population of dusty 
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starburst merging galaxies in our modelling, based on the near-IR number counts fitting 
and the background intensity. However, this value should be taken as an upper limit since 
the near-infrared light and stellar mass correlation will be biased in these starburst merging 
systems by the intense star formation activities, which can result in the near-infrared being 
dominated by young supergiants and the galaxies having a much smaller mass-to-light ratio 
than ellipticals. 

4.3. Redshift Distribution and Cosmic Background 

The redshift distribution of the three model components (starburst galaxies, "non- 
evolving" galaxies and dusty AGNs) with a detection limit SW^m > 0.083 mJy is presented 
in Fig. 9 (Left). A comparison with the recent photometric redshift estimation by Perez- 
Gonzdlez et al. (2005) shows that the starburst population in the modelling can only account 
for about 30% of the high-z tail (z > 2). At these redshifts, the strong PAH peak begins 
to move out of the MIPS 24/iin band and the strong mid-infrared continua from AGNs may 
play an increasing role in the detections. Therefore, we suspect that dusty AGNs, or the 
AGN population in a phase with luminosities still dominated by starbursts, may account for 
the additional objects in the high-z tail (Wang et al. 2003). We shall in the future properly 
handle such a population in expanded modelling. 

The right panel of Fig. 9 shows the calculated redshift distribution of the ISOCAM 15/im 
sources with the detection limit (> 0. 12 mJy). These luminous infrared sources cover a 
wide redshift range of 0.5 ~ 2, peaking at z ~ 1. Cowie et al. (2004) recently probed 
the evolution of the bright tail of the radio luminosity function up to z — 1.5 with redshift 
measurements of two deep 1.4 GHz fields. They found that the number density of the sources 
with ULIG radio power and no clear AGN signatures evolves as (1 + z) 7 up to about redshift 
z ~ 1. A comparison of the 1.4 GHz radio surveys with our model prediction on the redshift 
distribution of ULIGs is shown in the left panel of Fig. 10. We see that the number density 
of ULIGs increases dramatically up to redshift z ~ 1, but decreases afterwards. We know 
that the number density of starburst sources at z » 1 is still loosely constrained in the 
model, and recent detections of massive galaxies in the young Universe challenge most of the 
galaxy evolution theories. However, this part of discussion is out of the scope of our current 
work (Cimatti et al. 2004; Glazebrook et al. 2004). 

Finally, the infrared background from the model calculation is presented in Fig. 10 
(Right). Our predictions are consistent with the direct background measurements from 
COBE and other deep infrared surveys, as well as the upper limits from high energy TeV 7 
ray detections of nearby Blazars (Funk et al. 1998; Guy et al. 2000; Hauser & Dwek 2001). 
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5. Summary 

We have described a galaxy evolutionary scenario where starburst and AGN activities 
are triggered by galactic mergers/interactions. Assuming an infrared burst phase from gas 
rich mergers at high redshift, we have successfully interpreted the available observations 
ranging from near-infrared up to radio deep surveys, leaving the infrared background con- 
sistent with the upper limits from direct infrared background measurements by COBE and 
other infrared deep surveys, as well as recent TeV 7 ray detections of nearby Blazars. 

We further found from the current deep extragalactic surveys performed in the far- 
infrared with Spitzer, and at radio wavelengths with the VLA, that: 

1. ) LIGs with M82-like enhanced PAH emission features dominate the number counts at 
sub-mJy flux levels in Spitzer 24fim and ISOCAM 15/xm deep surveys; 

2. ) the 1.4 GHz source excess at sub-mJy levels and the strong evolution detected by ISO- 
CAM 15/im and Spitzer 24/im extragalactic deep surveys at flux S v < lmJy are due to 
the appearance of the same population, which is mostly a LIG phase triggered by galaxy 
interactions at z ~ 1; and 

3. ) considering near- and mid- infrared deep surveys and the background measurements, 
dusty starburst galaxies triggered by galactic interactions at z ~ 1 in the model account for 
a moderate fraction of the stellar mass in the Universe (30% to 50% as an upper limit). 

Although our modelling fits the full suite of available data well, it is not likely to 
be unique. To break the degeneracies in the model and make further progress requires 
both a better understanding of the infrared burst phase, and more accurate merger rate 
statistics. With the existing model constraints such as number counts, redshift distributions 
and background intensities, it is still very difficult to address such issues in detail and give 
an independent estimation of the model parameters. A real physical model including both 
proper chemical evolution and dust emission properties should be used in future work, for a 
better understanding of both the emission properties and dynamical processes of the dusty 
starburst merging galaxies across cosmic time. 
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Fig. 1. — Left: The SED of ULIGs selected for the modelling that gives a successful fit for 
the present deep surveys, especially in the far-infrared. The data points are photometric 
measurements of Arp220 (corrected to rest frame), collected by Elbaz et al. (2002); Right: 
A comparison of the SED of LIGs selected for the modelling with the photometric measure- 
ments of M82. The data points are collected by Silva et al. (1998). The meanings of the 
lines are denoted in the figure. 




Fig. 2. — Left: The model fitting of the differential number counts in Spitzer 24/im deep 
surveys. The data points are from Papovich et al. 2004, normalized to a Euclidean slope, 
which show the average counts from all the Spitzer fields and corrected for completeness; 
Right: The calculated contribution from LIGs, ULIGs and LIGs+ULIGs to the differential 
number counts at 24/xm. We see that LIGs dominate the Spitzer 24/im number counts at 
the sub-mJy flux level, while ULIGs contribute mostly to the bright part of the peak at 
0.3 m Jy < S u < lmJy. 
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Fig. 3. — Left: The model fitting of the differential number counts at 70/im. The blue 
asterisks (CDFS), blue empty diamonds (Marano field), and blue filled squares (Bootes field) 
are source counts in Spitzer deep surveys with no correction for incompleteness (Dole et al. 
2004). Red empty squares are IRAS 60/im counts from Lonsdale et al. (1990) converted 
to 70/im; Right: The differential number counts fitting at 160/xm by the model calculation. 
The blue stars (CDFS), blue empty diamonds (Marano field) are the source counts at 160/im 
with no correction for incompleteness (Dole et al. 2004). Blue filled squares correspond to 
ISO FIRBACK 170um number counts (Dole et al. 2001). 





Fig. 4. — Left: Model fitting of the differential number counts at ISOCAM 15 /im. The data 
points are from a variety of ISO deep surveys (Elbaz et al. 1999); Right: IRAS 60 \xm number 
counts fitting. The data points are from the IRAS Point Source Catalogue (1985) (PSC, 
black filled circles), Hacking et al. IRAS deep survey (HCH, red empty triangles), FSC (blue 
asterisks) from deep surveys by Moshir et al. (1992) and Saunders (1990). The kink seen at 
the bright end of the burst population number counts is due to the numerical inaccuracy of 
the logarithm binning during the Monte-Carlo simulation. 
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Fig. 5. — Left: The model fitting of the ELAIS normalized differential counts at 90 [im. 
The data are from recent 90/im final analysis by Heraudeau (2004) (empty squares) and 
Lockman Hole counts from Rodighiero (2003) (crosses), IRAS counts (stars) for galaxies in 
PSCz catalogue; Right: The differential number count fitting for FIRBACK 170 /zm deep 
survey by the model. The data points are from Dole et al. (2001). The kinks in the modelled 
number counts arise for the same reason as in Fig. 4. 
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Fig. 6. — Left: Model fitting of the integral number counts at 850 /im. The blue asterisks 
are from Blain et al. (1999), and the blue filled squares are from Barger et al. (1999); Right: 
The differential number count fitting of the 1.4 GHz radio deep survey by the modelling. 
The blue asterisks are the 1.4 GHz counts in the HDF-N from Richards (2000). The model 
result demonstrates that IR emitting starburst galaxies do not contribute significantly to the 
1.4 GHz counts at the bright end of the number counts diagram, but are important for the 
population of radio sources at flux densities < 0.5 — 0.8 m Jy. The 1.4 GHz number counts 
at the bright end (> 1 mJy) are primarily due to radio AGNs, which are not relevant to our 
discussion here. 
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Fig. 7. — Left: Model fitting of the cumulative number counts of infrared sources at 1.6/im 
brighter than a given threshold F v . The data points are from NICMOS observations covering 
1/8 — th of the Hubble Deep Field North, which reaches a flux level about three magnitude 
fainter than the deep surveys at K-band and ISOCAM 6.7/im (Thompson et al. 1999); Right: 
source count fitting by our modelling at 2.2/xm. Data are from McCracken et al. (2000). 
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Fig. 8. — Model fitting of the cumulative number counts of infrared sources at 6.7 fim with 
S u > 10 ji Jy within 16 arcmin 2 . Data are from Metcalfe et al. (2003, red empty triangles) 
and Sato et al.(2003, blue asterisks). 
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Fig. 9. — Left: The predicted redshift distribution by modelling of three major infrared 
populations, with a detection limit S24 Mm > 0.083 mJy. A comparison with the recent pho- 
tometric redshift estimation by Perez-Gonzdlez et al. (2005) is presented. Blue dot-dashed 
and red dotted thin lines are the measured upper and lower limits; Right: The predicted 
redshift distribution, within a detected limit > 0. 12 mJy at 15 //m. The histogram of the 
observed redshift distribution of ISOCAM 15/xm sources is taken from Franceschini et al. 
(2001). 




Fig. 10. — Left: The predicted number density distribution of ULIGs (dashed line) is com- 
pared with that of the starburst galaxies with ULIG radio power in the redshift bin 0.1 — 0.6, 
0.6 — 1.0 and 1.0 — 1.4 measured by Cowie et al. 2004 (blue filled squares). The blue aster- 
isks denote the local number density of ULIGs + LIGs and z = 1 — 3 radio/submm sources 
(> 6mJy) from Barger et al. 2000. The pink diamond shows the estimated value by Cowie 
et al. 2004 from the local star forming radio luminosity function; Right: Comparison of the 
calculated background spectrum with the measurements by independent groups in the all- 
sky COBE maps (e.g. Hauser et al. 1998), and the ultradeep optical extragalactic surveys 
by the HST in the HDF (black filled triangles, Madau & Pozzetti 2000). The three lower 
datapoints (red filled diamonds) in the far-IR are from re-analysis of the DIRBE data by 
Lagache et al.(1999), and the red fenced area is upper limits to the far-infrared background 
level from Fixsen et al. (1998). The two mid-IR points are the resolved fraction of the CIRB 
by the deep ISO surveys IGTES (empty red triangles, Elbaz et al. 2002b), while the dashed 
histograms are limits set by TeV cosmic opacity measurements. 



